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Abstract. The accurate modelling of astrophysical scenarios involving compact objects 
and magnetic fields, such as the collapse of rotating magnetized stars to black holes or 
the phenomenology of 7-ray bursts, requires the solution of the Einstein equations together 
with those of general-relativistic magnetohydrodynamics. We present a new numerical code 
developed to solve the full set of general-relativistic magnetohydrodynamics equations in 
a dynamical and arbitrary spacetime with high-resolution shock-capturing techniques on 
domains with adaptive mesh refinements. After a discussion of the equations solved and of 
the techniques employed, we present a series of testbeds carried out to validate the code and 
assess its accuracy. Such tests range from the solution of relativistic Riemann problems in 
flat spacetime, over to the stationary accretion onto a Schwarzschild black hole and up to the 
evolution of oscillating magnetized stars in equilibrium and constructed as consistent solutions 
of the coupled Einstein-Maxwell equations. 



PACS numbers: 04.25.Dm, 95.30.Qd, 04.40.Dg, 97.60Jd, 95.30.Sf 

1. Introduction 

Magnetic fields are ubiquitous in astrophysical objects and can play an important role, 
especially in those scenarios involving compact objects such as neutron stars and black holes. 
An accurate and consistent modelling of these scenarios, which are extreme both for the 
gravitational and the electromagnetic fields, cannot be done analytically and perturbative 
methods are also of limited validity. In the absence of symmetries, in fact, no dynamical 
and analytic solutions are known and it is only through the full solution of the equations 
of general-relativistic magnetohydrodynamics (GRMHD) that one can hope to improve our 
knowledge of these objects under realistic conditions. 

As in general-relativistic hydrodynamics, the work in this area of research has started, 
more than 30 years ago with the pioneering work of Wilson HI in lower spatial dimensions 
(see [2] for a review of the technical and scientific progress in relativistic hydrodynamics). 
Unlike general relativistic hydrodynamics, however, where both technical issues and scientific 
investigations have now reached an advanced stage of sophistication and accuracy, progress in 
GRMHD has yet to reach a comparable level of maturity. Indeed, it was only over the last few 
years that the slow but steady progress in GRMHD has seen a renewed burst of activity, with a 
number of groups developing a variety of numerical codes solving the equations of GRMHD 
under different approaches and approximations. This is partly due to the considerable added 
complexity of the set of equations to be solved in GRMHD and, partly, to the fact that only 
recently sufficient computational resources have become available to tackle this problem in 
two or three spatial dimensions and with sufficient resolution. 
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It does not come as a surprise, therefore, that most of the numerical codes developed in 
the last decade have been based on the same non-conservative formulation of the GRMHD 
equations introduced by Wilson, solving them on a fixed background to study accretion 
disks around black holes. In 0), for instance, the effects of a Kerr black hole on 
magnetohydrodynamical accretion have been studied, with particular attention being paid 
to the transfer of energy and angular momentum. Koide et al. 0, on the other hand, have 
developed a numerical code based on the artificial-viscosity approach proposed by Davis |6| 
to perform the first simulations of jet formation in General Relativity [7| and to study the 
possibility of extracting the rotational energy from a Kerr black hole OH]]. Furthermore, 
a distinct numerical code has been constructed by De Villiers and Hawley [ 1 1 using the 
formulation proposed in ATI [121 to carry out a series of studies on accretion flows around 
Kerr black holes lfl3l[T4l[T5ll. 

It is only rather recently that different groups have started to recast the system of 
GRMHD equations into a conservative form in order to benefit of the use of high-resolution 
shock-capturing schemes (HRSC) lfl6l [171 [181 [191 l20l . Such schemes, we recall, are essential 
for a correct representation of shocks, whose presence is expected in several astrophysical 
scenarios and in particular in those involving compact objects. Two mathematical results 
corroborate this view, with the first one stating that a stable scheme converges to a weak 
solution of the hydrodynamical equations EH , and with the second one showing that a non- 
conservative scheme will converge to the wrong weak solution in the presence of a shock |22|. 

All of these newly developed codes |16] [TTl |T8l [T9) that make use of HRSC methods 
have been so far applied to the study of accretion problems onto black holes, for which the 
self-gravity of the accreting material introduces very small corrections to the spacetime and 
fixed spacetime backgrounds can be used satisfactorily. 

Approaches alternative to that of constructing GRMHD codes have instead been based 
on the use of a modified Newtonian gravitational potential to mimic general relativistic 
effects without having to solve numerically Einstein equations (see [23] for an application 
to magnetorotational collapse of stellar cores) or on the use of different numerical methods, 
such as smoothed particle hydrodynamics and artificial viscosity, to study the merger of binary 
neutron star systems as a possible engine for short 7-ray bursts ll24l . Although the use of these 
approximations has made it possible to investigate this astrophysical scenario for the first 
time including details about the microphysics, it is clear that equally important corrections 
coming from the dynamical evolution of the spacetime need to be introduced when trying to 
model the phenomenology that is thought to be behind 7-ray bursts. As a first step in this 
direction, two codes were recently developed to solve the full set of GRMHD equations on 
a dynamical background l25ll26l . These codes, in particular, were used to perform the first 
study, in two spatial dimensions, of the collapse of magnetized differentially rotating neutron 
stars ll27l l28l l29l which are thought to be good candidates for short 7-ray bursts. 

Here, we present WhiskyMHD, a new three-dimensional numerical code in Cartesian 
coordinates developed to solve the full set of GRMHD equation on a dynamical background. 
The code is based on the use of high-resolution shock-capturing techniques on domains with 
adaptive mesh refinements, following an approach already implemented with success in the 
general-relativistic hydrodynamics code Whisky l30l . and which has been used in the study 
of several astrophysical scenario with particular attention to gravitational-wave emission from 
compact objects. 

The paper is organized as follows: in Sect.|2]we recall the equations of GRMHD and the 
form they assume when recast in a conservative form, while in Sect. [3] we discuss in detail the 
numerical methods adopted for their solution. Section|4]is dedicated to the series of testbeds 
the code has passed both in special and in general relativistic conditions. Finally, Sect. goffers 
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a summary of the results and an overview on our future projects. 

Throughout the paper we use a spacelike signature (—,+,+,+) and a system of units in 
which c = G = Mq = 1. Greek indices are taken to run from to 3, Latin indices from 1 to 
3 and we adopt the standard convention for the summation over repeated indices. Finally we 
indicate 3-vectors with an arrow and use bold letters to denote 4-vectors and tensors. 



2. Formulation of the equations 

We adopt the usual 3 -dimensional foliation of the space time so that the line element reads 

ds 2 = -(a 2 - f3 l Pi)dt 2 + 2/Wdi + j lj dx i dx j , (1) 

where (3 % is the shift vector, a is the lapse function and Yy are the spatial components of the 
four-metric g M „. 

As its predecessor Whisky 11301 . the WhiskyMHD code benefits of the Cactus 
computational toolkit RD which provides an infrastructure for the parallelization and the 
I/O of the code, together with several methods for the solution of the Einstein equations. As 
a result, at each timestep our new code solves the MHD equations while Cactus provides the 
evolution of the metric quantities. The evolution of the field components is done using the 
NOK formulation [32, 33, 34 1 and details about its numerical implementation can be found 



Here too we make use of the so-called "Valencia formulation" 11381 [39ll which was 
originally developed as a 3 + 1 conservative Eulerian formulation of the general relativistic 
hydrodynamic equations, but which has been recently extended to the case of GRMHD [ 18 1. 
Following [18 1 we define the Eulerian observer as the one moving with four velocity n 
perpendicular to the hypersurfaces of constant t at each event in the spacetime. This observer 
measures the following three-velocity of the fluid 



-u^n^ W a 



(2) 



where = + n^riy is the projector orthogonal to n, u is the four-velocity of the 
fluid and —u^n^ = au° = W is the Lorentz factor which satisfies the usual relation 
W — 1/VT — v 2 , where v 2 = YyuV. The covariant components of the three-velocity 
are simply given by Vi — Ui/W . 

2.1. Maxwell equations 

The electromagnetic field is completely described by the Faraday electromagnetic tensor field 
F^ v obeying Maxwell equations (cfr ll40l ) 

V„ T» v = , (3) 

\7 U F^ = 4tt J* , (4) 

where V is the covariant derivative with respect to the four-metric g, J' is the charge current 
four-vector and *F is the dual of the electromagnetic tensor defined as 

= ^ X5 F X5 , (5) 

rj^ uXS being the Levi-Civita pseudo-tensor. A generic observer with four-velocity U will 
measure a magnetic field B and an electric field E given by 

E a = F al3 Uf3 , (6) 

B a = *F afs U , (7) 
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and the charge current four-vector J' can be in general expressed as 

J* = qu» + o-F^Uu , (8) 

where q is the proper charge density and a is the electric conductivity. 

Hereafter we will assume that our fluid is a perfect conductor (ideal MHD limit) and so 
that o — > oo and F^ v u v — (i.e. the electric field measured by the comoving observer is 
zero) in order to keep the current finite. In this limit, the electromagnetic tensor and its dual 
can be written exclusively in terms of the magnetic field b measured in the comoving frame 

puc = v an»* baU ^ ; * Ffl v = h n u v _ h v ul x (9) 

with the Maxwell equations taking the simple form 

V„ *F» V = -^=d v ( (&"«" - = , (10) 

In order to express these equations in terms of quantities measured by an Eulerian observer, 
we need to compute the relation between the magnetic field measured by the comoving and by 
the Eulerian observers, respectively b and B. To do that we introduce the projection operator 
= + u^Uv orthogonal to u. If we apply this operator to the definition of the magnetic 
field B measured by an Eulerian observer, we can easily derive the following relations 

a ' W M VF 2 ' 

where B 2 = B l Bi. The time component of equations ( fTOb provides the divergence-free 
condition 

d t B % = , (12) 

where B l = ^/jB l and 7 is the determinant of 7^ . The spatial components of equations ( fTol i. 
on the other hand, yield the induction equations for the evolution of the magnetic field 

d t (B l ) = dj(v l B j — v^B i ) , (13) 

where v l = av l — (3 l . 

2.2. Conservation equations 

The evolution equations for the rest-mass density p, the specific internal energy e and for 
the three-velocity v l can be computed, as done in relativistic hydrodynamics, from the 
conservation of the baryon number 

V v {pu v ) = , (14) 

from the conservation of the energy-momentum 

V„T^ = , (15) 

and from an equation of state (EOS) relating the gas pressure p to the rest-mass density p 
and to the specific internal energy e. We also assume that the fluid is perfect so that the total 
energy momentum tensor, including the contribution from the magnetic field, is given by 

= {ph + b 2 ) u»u v + (p + ^ g^ v - b% v , (16) 

where h^l + e + p/pis the specific relativistic enthalpy. 
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Following lfl8ll and in order to make use of HRSC methods, we rewrite equations (TBI . 
< n~5b and ( fT3l ) in the following conservative form 

[a t (V7F°) + diW^gF)] = S , (17) 

where F° = (D, Sj, t, B k ) T is the vector of the conserved variables measured by the Eulerian 
observer, F l are the fluxes 

/ Dtf/a \ 

SjV l /a +(p + b 2 /2) 8) - bjB l /W 
TV 1 /a +(p + b 2 /2) v i - ab^/W 

B k v i /a - B i i k /a J 



\ 

and S are the source terms 



(18) 



where 



\ 

TV (dtfvj - T 5 vll g Sj ) 
aiT^d^a-T^T ^) 

Qk 



D = pW , 

Sj = (ph + b 2 )W 2 v j - OLb% 3 , 

t ^(ph + b 2 )W 2 ~(p+^)-a 2 (b ) 2 



(19) 



D 



(20) 
(21) 

(22) 



and k — (0, 0, 0) T . While ready to make use of arbitrary EOS, the ones presently 
implemented in the code have a rather simple form and are given by either the polytropic 
EOS 



p = Kp L 



P 



p(T-l)> 
or by the ideal-fluid EOS 

p=(T-l)pe, 

where K is the polytropic constant and T is the adiabatic exponent. In the case of the 
polytropic EOS ( l23l ). V = 1 + 1/N, where N is the polytropic index. 



(23) 
(24) 

(25) 



3. Numerical methods 



As in the Whisky code, the evolution equations are integrated in time using the method of 
lines [41 1, which reduces the partial differential equations (fTTT i to a set of ordinary differential 
equations that can be evolved using standard numerical methods, such as Runge-Kutta or the 
iterative Cranck-Nicholson schemes ll42ll43l . 
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3.1. Approximate Riemann solver 

As mentioned in the Introduction, WhiskyMHD makes use of HRSC schemes based on 
the use of Riemann solvers to compute the fluxes between the numerical cells. More 
specifically, we have implemented the Harten-Lax-van Leer-Einfeldt (HLLE) approximate 
Riemann solver |44| which is simply based on the calculation of the eigenvalues of eqs. ( fTTI ) 
and it does not require a basis of eigenvectors. In the HLLE formulation the flux at the 
interface between two numerical cells is therefore computed as 



where PJf and are computed from the values of the primitive variables reconstructed at 
the right and left side of the interface P r and Pi, respectively. The coefficients c max = 
max(0, c+, r , c+,;)> c m ; n = — min(0, c_ ir , c_.z) and c± jr , c±.i are instead the maximum left- 
and right-going wave speeds computed from P r and Pi. In our implementation P r and 
Pi are computed using a second order TVD slope limited method which can be used with 
different limiters such as minmod, Van Leer and MC |41 1. 

An alternative to the use of an approximate Riemann solver could have been the use 
of the exact Riemann solver recently developed in GRMHD. We recall that in relativistic 
hydrodynamics the exact solution is found after expressing all of the quantities behind each 
wave as functions of the value of the unknown gas pressure p at the contact discontinuity. 
In this way, the problem is reduced to the search for the value of the pressure that satisfies 
the jump conditions at the contact discontinuity (see l45l l46l l47l |48 1 for the details). The 
procedure for the exact solution of the Riemann problem in relativistic MHD is based instead 
on the use of an hybrid approach that makes use of different set of unknowns depending on the 
wave. In the case of fast-magnetosonic waves, both shocks or rarefactions, all the variables 
behind the waves are rewritten as functions of the total pressure, i.e. p + b 2 /2, while behind 
slow magnetosonic shocks or rarefactions the components of the magnetic field tangential to 
the discontinuity are used to compute all the other variables. The use of this strategy was 
essential in order to reduce the problem to the solution of a closed system of equations that 
can be solved with standard numerical routines such as Newton-Raphson schemes (see [49 ] 
for the details). The numerical code computing the exact solution is freely available from the 
authors upon request and it is now becoming a standard tool in the testing of both special and 
general relativistic MHD codes. 

While the use of an exact Riemann solver could provide the solution of the discontinuous 
flow at each cell interface with arbitrary accuracy, the exact solver presented in |49| is still 
computationally too expensive to be used in ordinary multidimensional codes and we have 
found the HLLE algorithm to be sufficiently accurate for the resolution used in our tests. 

3.1.1. Calculation of the eigenvalues An important difference with respect to relativistic- 
hydrodynamic codes is that in GRMHD the calculation of the eigenvalues required by 
the HLLE solver is made more complicated by the solution of a quartic equation. The 
characteristic structure of GRMHD equations is analyzed in detail in Boll and we simply 
report here the expressions for the calculation of the seven wave speeds associated with the 
entropic, the Alfven, the fast and slow magnetosonic waves. 

More specifically, the characteristic speed A of the entropic waves is simply given by 
v = av l — f3 l , while the values for the left- and right-going Alfven waves are 



F l = 




(26) 




b l ± ^{ph + b 2 )u l 



(27) 



6° ± y/(ph + b 2 )u° 
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Similarly, the four speeds that are associated with the fast and slow magnetosonic waves, 
and that are required in the calculation of the fluxes, can be obtained by the solution of the 
following quartic equation in each direction i for the unknown A 

Ph(^-l)a 4 - (ph + ^]a 2 G + B 2 G = , (28) 

where 



ir 



X + av'-P) , (29) 
a " 

B = V - b°\ , (30) 

G= \ [-(A + /? ! ) 2 +aV l ] , (3D 

and c s is the sound speed (Note that the convention on repeated indices should not be used 
for the last term in the expression for G, i.e. 7"). In the degenerate case in which B % = 0, 
eq. (f28b can be reduced to a simple quadratic equation that is solved analytically. In the more 
general case, however, eq. (|28| | cannot be reduced to the product of two quadratic equations 
as in Newtonian MHD and different methods are implemented in the code in order to solve 
it. The first one simply makes use of the analytic expression for a quartic equation |50|, while 
the other two search the solution numerically either through an eigenvalue method or through 
a Newton-Raphson iteration BP . The latter has shown to be the most accurate and robust and 
it is the one used by default. 

We have also implemented an approximate method for the calculation of the eigenvalues 
associated with the fast magnetosonic waves (which are the only two roots needed by HLLE) 
which was introduced in [52 1 and which reduces the original quartic to a quadratic equation, 
that can be solved analytically, by imposing B % = and B^Vj — in equation (1281 l, The 
values computed in this way differ by less than 1% with respect to the exact values and we 
have used them in those situations in which the solution of the quartic can be complicated by 
the presence of degeneracies or when two of the roots are very close to each other. 



3.2. Constrained-Transport Scheme 

Although an exact solution of eqs. dT3l would guarantee that the constraint condition (V2\ is 
also satisfied identically and all times, any numerical solution of the induction equations dT3l l 
will inevitably produce a violation of the divergence-free condition which, in turn, may lead 
to unphysical results or even to the development of instabilities [53|. To avoid this problem 
several numerical methods were developed in the past starting from the so-called "staggered 
mesh magnetic field transport algorithm" first proposed by Yee [54] and then implemented in 
an artificial-viscosity scheme with the name of "constrained-transport" scheme (CT) by Evans 
&Hawley (55]. 

A modified version of the CT scheme ll55l which is based on the use of the fluxes 
computed with HRSC methods has been proposed by Balsara & Spicer [56 1 ("flux-CT") and 
is the one implemented in our code because of its simplicity and computational efficiency. An 
interested reader is referred to [57] for other possible methods to enforce the divergence-free 
condition with HRSC schemes. 

We recall that the flux-CT method is based on the relation that exists in ideal MHD 
between the fluxes of the magnetic field B and the value of the electric field E = —v x B. In 
particular, if we define F 4 = a^FfF' 1 then the following relations hold 

E x = F z (§ yS j = -F v (b z ^j , (32) 
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EP> = -F z = F x (B z ^j , (33) 

E z = F y (B x ^ = -F x (bA , (34) 

where F l (b^ = v % B^ — xPB % . The induction equation ( fT3l l can then be written as 

d t B + V x E = . (35) 

Taking the surface integral of ( 1351 ) across a surface £ between two numerical cells, Stokes' 
theorem yields 



d t / B • dS + / E ■ I = , (36) 

where Z is the unit vector parallel to the surface boundary <9£. Considering for simplicity the 
^-direction, with the surface £ as located at (i + k) and the integers (i, j, k) denoting the 
cell centers on our discrete grid (see figure[T|i, we can define 

BLl , fc = f B-dt, (37) 



i+Tid,k AyAz „ . 
and use the finite-difference expression of eq. d36b to obtain 

dtBf.l I =- (e V , ! .. 1 -E V , ! . Ull ) /AZ- (E Z 1 ! . -Sf, ! - 1 ,) /Ay, 

(38) 

where the values of the electric field on the edges of the surface are simply computed taking 
the arithmetic mean of the fluxes across the surfaces that have that edge in common [cf (f32l> — 

mi e.g. 

E l+\ ,j,fe+± = 4 ( F i+±,j,k + ^+ij,fc+i ~ ^j,fc+i _ ^-ij.fc+i) ^ (39) 

where the fluxes F l are those computed with the approximate Riemann solver. 

Since we are using HRSC methods, all the quantities are located at cells centers but in 
equation (l38l we are effectively evolving the magnetic field at the surfaces between the cells. 
The relation between these two different values of the magnetic field is given by a simple 
average 



^ = 5(^ + ^-4,0 • ^ 

• -) ■ (42) 

To demonstrate that this method guarantees that V • B = and will not grow in time, 
we can integrate over the volume of a numerical cell this constraint and then using the Gauss 
theorem we obtain 

6 r - 

V ■ BdV = 22 B-dS, (43) 

AV i=1 Js, 

where the sum is taken over all the six faces S, that surround the cell. Taking now the time 
derivative of this expression and using eq. (f36t we obtain 

d t I V ■ BdV = -Y\j E-T, (44) 
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Figure 1. Schematic diagram of the data needed for the CT scheme. The evolution of 
j k ' s determined by the values of the electric field E at the edges of the surface 
S located at (i + 1/2, j, ft), 



and the sum on the right-hand side gives exactly zero since the value of E ■ I for the common 
edge of two adjacent faces has a different sign. 

3.3. Primitive-variables recovery 

Because the fluxes F 2 depend on the primitive variables P and not on the evolved conservative 
variables F°, the values for the primitive variables need to be recovered after each timestep 
and at each gridpoint. With the exception of the magnetic field variables B l , the complexity 
of the system of equations to be solved prevents from an analytic solution relating the 
primitive to the conservative variables through simple algebraic relations and thus the system 
of equations d20ii-(l22ii needs to be solved numerically. Several methods are available for 
this, the most obvious (and expensive) one consisting of solving the full set of 5 equations 
given by the expressions for (D, Sj, t) in the 5 unknowns (p, v l , e); we refer to this as to 
the 5D method. Alternatively, and under certain conditions, it is possible to reduce the set of 
equations to be solved to a couple of nonlinear equations (2D method) or even to a single one 
(ID method). We review them briefly in the following Sections but a more detailed discussion 
can be found in |58|. 

3.3.1. 2D method The following procedure is the same used in iTTSll and it is an extension to 
full General Relativity of the method developed in |59| in special relativity. The idea is to take 
the modulus S 2 = 5 J Sj of the momentum instead of the expression for its three components 
reducing the total number of equations that one has to solve. Using the relations (fTTl i it is 
possible to write S 2 as 

S 2 = (Z + B 2 ) 2 ^-——^- - (2Z + B 2 ) { -^f- , (45) 



WhiskyMHD: a new numerical code for general relativistic magnetohydrodynamics 10 



where Z = phW 2 . It is also possible to rewrite the equation for the total energy in a similar 
way 

Using then the definition of D = pW, eqs. d45l l and (l46*l l form a closed system for the two 
unknowns p and W, assuming the function h = h(p,p) is provided. When using a polytropic 
EOS [i.e. eq. (1231) 1. the integration of the total energy equation is not necessary (the energy 
density can be computed algebraically from other quantities) and the system reduces to the 
numerical solution of the equation d45t . Once the roots for W, p and p = D/W are found, it 
is possible to compute the values of vt using the definition of the momentum Si 

_ BjjB^+SjZ 



Z(B 2 + Z) 



(47) 



3.3.2. ID method The basic idea of this method is to consider also the gas pressure p as a 
function of W reducing the total number of equations that must be solved numerically. When 
using an ideal-fluid [i.e. eq. d25lll. Z can in fact be rewritten as 

Z = DW + -^—^p{W)W 2 . (48) 

Using equation d48l ) it is possible to rewrite ( |46T > as a cubic equation for p(W) which admits 
only one physical solution. So at the end we need only to solve equation (|43T > for the only 
unknown W . Having obtained W, we can then compute p = p( W) and the other quantities 
in the same way as done in the 2D method. 



3.4. Atmosphere treatment 

As already done in the Whisky code and in other full GRMHD codes ll26l l25ll we avoid 
the presence of vacuum regions in our domain by imposing a floor value to the rest-mass 
density. This is necessary because the routines that recover the primitive variables from the 
conservative ones may fail to find a physical solution if the rest-mass density is too small. 
The floor value used for the tests reported here is p atm = 10~ 7 x max(po), with po being the 
value of the rest-mass density at t = 0, but a floor which is two orders of magnitude smaller 
works equally well. In practice, for all of the numerical cells at which p < p a t m , we simply 
set p = p a tm, v-' — 0, and do not modify the magnetic field. This is different from what done 
in other codes (e.g., (26] |25]), which set to zero the initial value of the magnetic field in the 
low density regions. 



3.5. Excision 

Many interesting astrophysical scenarios involve the presence of black holes and so of regions 
of spacetime where singularities are present. These regions are causally disconnected from 
the rest of the physical domain and the values of the fields inside should not affect the zone 
outside the event horizon. This is not true in numerical codes where it can happen that some 
information from inside the event horizon is used to compute the values of the variables 
outside. In order to avoid this, excision algorithms were developed in general relativistic 
hydrodynamics and they are based on the use of some kind of boundary condition applied 
to the boundary between the excised zone, where the equations are no more solved, and the 
domain outside. As already done in the Whisky code we apply a zeroth-order extrapolation 
to all the variables at the boundary, i.e. a simple copy of the MHD variables across the 
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Table 1. Initial conditions for the Riemann problems used to test the code. 



Test type 


State 


P 


P 


v x 




v z 


B x 


B y 


B z 


Balsara Test 1 




















(T = 2) 


left 


1.000 


1.0 


0.0 


0.0 


0.0 


0.5 


1.0 


0.0 




right 


0.125 


0.1 


0.0 


0.0 


0.0 


0.5 


-1.0 


0.0 


Balsara Test 2 




















(r = 5/3) 


left 


1.0 


30.0 


0.0 


0.0 


0.0 


5.0 


6.0 


6.0 




right 


1.0 


1.0 


0.0 


0.0 


0.0 


5.0 


0.7 


0.7 


Balsara Test 3 




















(r = 5/3) 


left 


1.0 


1000.0 


0.0 


0.0 


0.0 


10.0 


7.0 


7.0 




right 


1.0 


0.1 


0.0 


0.0 


0.0 


10.0 


0.7 


0.7 


Balsara Test 4 




















(T = 5/3) 


left 


1.0 


0.1 


0.999 


0.0 


0.0 


10.0 


7.0 


7.0 
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excision boundary. A different method, based on the use of a linear extrapolation, has instead 
been implemented in ||251 and although it can lead to improved accuracy for smooth flows 
(and especially when the MHD variables change rapidly near the excision boundary), it also 
leads to significantly incorrect results when shocks are present (see Sect. 14.21 ). Great care 
must therefore be paid at the properties of the flow near the excision boundary and the code 
presently includes both algorithms. 

It is important also to note that other methods, not based on excision techniques, are 
being developed to improve the stability of numerical codes when black hole are present in 
the domain. One of these approaches is based on the use of a Kreiss-Oliger dissipation for the 
field variables inside the apparent horizon 11601 and it can be straightforwardly extended also 
to the MHD case. 

3.6. Mesh Refinement 

The developments made in Wh i s ky for handling non-uniform grids have been extended also 
to Whi skyMHD which can therefore use a "box-in-box" mesh refinement strategy |61 1. This 
allows to reduces the influence of inaccurate boundary conditions at the outer boundaries and 
for the wave-zone to be included in the computational domain. In practice, we have adopted 
a Berger-Oliger prescription for the refinement of meshes on different levels |62| and used 
the numerical infrastructure described in lETTl . i.e., the Carpet mesh refinement driver for 
Cactus (see [63] for details). In addition to this, a simplified form of adaptivity is also 
in place and which allows for new refined levels to be added at predefined positions during 
the evolution or for refinement boxes to be moved across the domain to follow, for instance, 
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regions where higher resolution is needed. 

4. Tests 

Code-testing represents an important aspect of the development of any newly developed 
and multidimensional code because it validates that all of the algorithms are implemented 
correctly and represent a faithful and discretized representation of the continuum equations 
they are solving. In what follows we report the results for a series of testbeds ranging from 
the solution of relativistic Riemann problems in flat spacetime, over to the stationary accretion 
onto a Schwarzschild black hole and up to the evolution of isolated and oscillating magnetized 
stars. We note that in the tests involving a poly tropic EOS the recovery of the primitive 
variables has been made using the "2D-method" (see Sect. I3.3.T1 ). while a "ID-method" has 
been used when adopting an ideal-fluid EOS (see Sect. l3.3.2l >. 

4.1. Riemann problems 

As customary in the testing of hydrodynamics and magnetohydrodynamics codes we have 
first validated WhiskyMHD against a set of Riemann problems in a Minkowski spacetime 
following the series of initial conditions proposed by Balsara (64). All these tests were run on 
a grid of unit length with 1600 grid points with the initial discontinuity located at the center 
of the grid. An ideal equation of state with T = 5/3 was used with the exception of the first 
test with r = 2 and the initial conditions for all the tests are reported in Table [T] 

In all of the tests presented here the numerical solution for the different MHD variables 
has been compared with the exact one computed with the exact Riemann solver discussed 
in (49). This represents an important difference with what done in the past by similar codes 
as it allows, for the first time, for a quantitative assessment of the code's ability to evolve 
correctly all the different waves that can form in relativistic MHD. In figures[2]and[3]the exact 
solution is represented with a solid line, while the numerical one with different symbols. 

In figure [2] in particular, we show the comparison between the numerical and the exact 
solution at t = 0.4 for several MHD variables as computed for the relativistic analogue of the 
classical Brio-Wu shock tube problem 1 65 , 66 ] . The initial discontinuity develops a left-going 
fast rarefaction, a left-going slow compound wave, a contact discontinuity, a right-going slow 
shock and fast rarefaction. Note that besides presenting the solution with a = 1 and (3 X = 0, 
we have exploited the freedom in choosing these gauges and validated the code also for less 
trivial values of the lapse and shift |[T8ll . More specifically, shown with different symbols 
in figure[2]are the numerical solutions with a = 2 at time t = 0.2 and with (3 X = 0.4 after 
the latter has been shifted in space by f3 x t. Clearly, all the symbols overlap extremely well, 
coinciding with the exact solution also in the presence of strong discontinuities. Note that 
only 160 of the 1600 data points used in the simulation are shown and that the difference 
between the numerical solution and the exact one at the compound wave is due to the fact 
that, by construction, our exact solver assumes compound waves never form. We have indeed 
adopted the same standpoint of Ryu and Jones [ 67 1 in the development of their exact Riemann 
solver in nonrelativistic magnetohydrodynamics. We also remark that it is not yet clear 
whether compound waves have to be considered acceptable physical solutions of the ideal 
MHD equations and a debate on this is still ongoing (see, for instance, [68, 69] ESI EH)- 

Similarly, shown in figure [3] are the comparisons between the numerical solution for the 
rest-mass density and the T/-component of the magnetic field for the tests number 2 (first row), 
3 (second row), 4 (third row) and 5 (fourth row) of Balsara. The first 3 are computed at a 
time t — 0.4 while the test number 5 is computed at t = 0.55. Clearly, our code is able 
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Figure 2. Numerical solution of the test number 1 of Balsara with different values for the lapse 
a and the shift f3 x . The solid line represents the exact solution, the crosses the numerical one 
at time t = 0.4, the open triangle at time t = 0.2 but with a = 2 and the open squares at 
t = 0.4 but with f3 x = 0.4, in this last case the solution is shifted on the x-axis by the amount 
j3 x t. Note that only 160 of the 1600 data points used in the numerical solution are shown. 
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Figure 3. Numerical solution of the tests number 2 (first row), 3 (second row), 4 (third row) 
and 5 (fourth row) of Balsara. The first 3 are computed at a time ( = 0.4 while the test number 
5 is computed at t = 0.55. The solid line represent the exact solution while the open squares 
the numerical one. Note that only 160 of the 1600 data points used in the numerical solution 
are shown. 
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Figure 4. Numerical solution of the test number 2 of Balsara at a time t = . 4 with an excision 
boundary (dashed vertical line) located at x = 0.25; the region at the right of this boundary is 
not evolved. In the two left panels a zeroth-order extrapolation, i.e. a simple copy, was used, 
while in the two right panels the values of the different variables at the excision boundary 
were obtained with a linear extrapolation. The solid line represents the exact solution while 
the open squares the numerical one. The solution is composed of two left-going fast and slow 
rarefactions, of a contact discontinuity and of two right-going fast and slow shocks. Only 160 
of the 1600 data points used in the numerical solution are shown. 



to resolve all the different waves present in MHD, showing a very good agreement with the 
exact solution. Other Riemann problems have been carried out in different directions (either 
along coordinate axes or along main diagonals) and they all provide the same level of accuracy 
discussed in figures|2]and[3] 

4.2. Excision tests on aflat background 

We next show the code's ability to accurately evolve shocks also when an excised region is 
present in the domain. To this scope, we have used the test number 2 of Balsara excising 
the region x % £ [0.25,0.5] and using the zeroth-order extrapolation scheme. In this case 
the fast and slow shocks moving to the right go inside the excised region and the solution 
outside is not affected. This is shown in the left panels of figure [4] which report with small 
squares the numerical solution for p and B y of the test number 2 of Balsara with an excision 
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boundary (dashed vertical line) located at x = 0.25. The data refers to time t = 0.4, when the 
right-going waves have already gone through the excision boundary as indicated by the exact 
solution (continuous line). 

As a comparison, and to underline its incorrectness in the case of non-smooth flows, 
we have also considered boundary conditions involving a linear extrapolation of the MHD 
variables across of the excision boundary as suggested in ll25ll . This is illustrated in the right 
panels of figure 2] which are the same as the left ones but for the different boundary condition 
at the excision boundary. Clearly, in this case the solution outside the excision region is badly 
affected and a left-going wave is produced which rapidly spoils the solution. Because this 
happens only when the discontinuity crosses the excision boundary, it is clear that a linear 
extrapolation is not adequate in this case as it provides an incorrect information on the causal 
structure of the flow near the boundary. As we will discuss in the following Section, however, 
a linear extrapolation remains a good, and sometimes preferable, choice in the case of smooth 
flows. 

4.3. Magnetized spherical accretion 

This second test proves the ability of the code to evolve accurately stationary accretion 
solution in a curved but fixed spacetime. In particular, we consider the spherical accretion of 
a perfect fluid with a radial magnetic field onto a Schwarzschild black hole (this is sometimes 
referred to as a relativistic Bondi flow). The solution to this problem is already known for the 
unmagnetized case, but it is simple to show that its form is not affected if a radial magnetic 
field is added [10]. The initial setup for this test is the same used in iflOl [161 l25l [18 1 and 
consists of a perfect fluid obeying a poly tropic EOS with T = 4/3. The critical radius of 
the solution is located at r c = 8M and the rest-mass density at r c is p c — 6.25 x 10~ 2 . 
These parameters are sufficient to provide the full description of the accretion onto a solar 
mass Schwarzschild black hole as described in [72]. We solve the problem on a Cartesian 
grid going from x l = to a; 1 = 11M. 

To avoid problems at the horizon, located at r = 2M, the metric is written in terms of 
ingoing Eddington-Finkelstein coordinates. The excision boundary has the shape of a cubical 
box of length M so that the domain [0, M\ x [0, M] x [0, M] is excluded from the evolution. 
Furthermore, as a boundary condition across the excised cube we have considered both a 
zero-th and a first-order extrapolation, finding the latter to yield sligthly more accurate results 
(e.g., the overall error is smaller of w 3% for a test case with b 2 / p — 25, where b 2 j p is the 
dimensionless magnetic field strength as measured at r — 2M). As discussed earlier, this is 
indeed to be expected for smooth flows as the ones considered here for the relativistic Bondi 
flow. 

We measure the order of accuracy of the code by using the Li-norm of the relative error 
on the rest-mass density 

= Ei.j.fc \Pi,j,k - Pcxactfoi, Vj,Z k )\ 

We plot this quantity in the left panel of figure [5] as a function of the magnetic-field strength 
and as computed at time t = 100AT for two different resolutions of 100 3 and 150 3 gridpoints, 
respectively. In addition, the error from the high-resolution simulation is multiplied by 1.5 2 
so that the two curves should overlap if the code were second-order convergent. Clearly, the 
code does not show the expected convergence rate but for relatively weak magnetizations, 
i.e. b 2 /p < 4 (we recall that these corresponds nevertheless to rather large magnetic fields of 
ps 10 19 G). This behaviour is indeed similar to what found by Duez et al. [25 1 and has a rather 
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Figure 5. Left panel: L\ norm of the relative error in the rest-mass density for the magnetized 
spherical accretion test, shown for different values of the magnetic field. Results from 100 3 
and 150 3 runs are compared at time t = 100M, with the high-resolution curve being 
multiplied by 1.5 2 so that with a second-order convergence the two lines would overlap. Right 
panel: Relative error computed at point x = AM, y = 0, z = 0; also in this case the high- 
resolution line is scaled to produce an overlap in the case of second-order convergence. In both 
cases the insets offer a magnification for small values of the magnetic field. 



simple explanation. It is due to the rather large error in gridpoints near the excision boundary, 
i.e., for x l € [M, 2M), which spoil the overall behaviour of the L\ norm (admittedly not a 
good measure of the convergence for a solution which is so rapidly varying near the excision 
boundary). To clarify this, we show in the right panel of figure[5]the same as in the left panel 
but for the relative error computed at a single gridpoint, i.e., at x = 4M, y = 0, z = 0. Clearly 
the convergence is much closer to second-order in this case (the precise order being w 1.8) 
and for much larger range in magnetizations. 

A closer look at the behaviour of the relative error is offered in figure [6] where it 
is shown as measured along the x-direction for a magnetization of b 2 / p = 25 (i.e., with 
Pmag/Pgas = 97) and at time t = lOOAf. Also in this case, the high-resolution relative error 
is multiplied by 1.5 2 so that the two lines overlap if second-order convergent. Clearly this 
does not happen but also only for a small number of gridpoints near the excision boundary 
located at x = M. 

As a final remark, we point out that the simulations of spherical-accretion flows 
performed here span a range of magnetizations well beyond what considered in the past with 
codes using Cartesian coordinates (the results reported in l25ll . for instance, were limited to 
b 2 /p < 30). Indeed, no sign of instability has been found and only a moderate loss of accuracy 
has been measured for magnetic fields as large as b 2 / p < 160. 

4.4. Evolution of a stable magnetized neutron star 

As a final test validating the code in a fully dynamical evolution of both the MHD variables 
and of the spacetime, we now consider the evolution of a stable magnetized neutron star. 
Although this initial data represents a stationary solution, small oscillations can be triggered 
by the small but nonzero truncation error. Such oscillations are sometimes considered a 
nuisance and even suppressed through the introduction of artificial-viscosity terms. On the 
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Figure 6. Relative error of the rest-mass density for the magnetized spherical-accretion flow 
with b 2 /p = 25 (i.e., with Pmag Jpgas — 97) along the a>axis at t — lOOAf. The high- 
resolution error is multiplied by 1.5 2 so that the two lines would overlap with a second-order 
convergence. Clearly, this does not happen near the excision boundary located at x = M. 



contrary, since they represent the consistent response of the star to small perturbations, they 
should considered as extremely useful. The eigenfunctions and eigenfrequencies of these 
oscillations, in fact, can serve both as a test of the code, when compared with the expectations 
coming from perturbation theory (see Appendix B of 1 37 1 for a representative example), or 
to extract information on the properties of the star, when considering regimes which are not 
yet accessible to perturbative studies (e.g., in the case of nonlinear oscillations or very rapidly 
rotating stars). 

Two options are possible for the construction of the initial data. A first and simpler 
one was employed extensively in ll27l l28l |29l and consists of considering a background 
purely-hydrodynamical solution in stable dynamical equilibrium and of "adding" a poloidal 
magnetic field in terms of a purely toroidal vector potential. Besides being essentially 
arbitrary, the vector potential is chosen to be proportional to the pressure so as to lead to 
a magnetic field entirely confined within the star. While straightforward, this approach does 
not construct a magnetized stellar model which is consistent solution of the Einstein equations 
and thus inevitably introduces violations of the Hamiltonian and momentum constraints. Such 
violations, however, are in general negligible for reasonably small magnetic fields. 

A second option, and the one employed here, consists of computing the initial conditions 
as a consistent and accurate solution of the Einstein equations for a stationary, axisymmetric 
and magnetized star. We have done this by using the spectral code developed by Bocquet 
et al. [73], which solves the full set of Einstein and Maxwell equations to high precision. 
Assuming an axisymmetric model with a poloidal magnetic field having the dipole moment 
aligned along the rotation axis, the code is used to build initial configurations of uniformly 
rotating magnetized neutron stars with different angular velocities and magnetic field 
strengths. 
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Magnetic field 




For simplicity we here consider a nonrotating magnetized neutron star with mass M = 
1.3-Mq endowed with a poloidal magnetic field with magnetic dipole along the z-axis and a 
central magnetic field B c — 2.4 x 10 14 G, corresponding to (3 = Pmag/.Pgas = 10~ 6 (this 
(3, which should not be confused with the shift vector (3 l , is always monotonically decreasing 
inside the initial equilibrium model, and is much larger in the atmosphere, where it reaches 
values ~ 10 6 ). A polytropic equation of state with T = 2 and K = 372 was used both for 
the calculation of the initial model and during the evolution. A representative image of the 
initial model is presented in figure [7] which shows the magnetic field lines together with the 
stellar surface (thick solid line). Note that although the star is nonrotating, the presence of a 
magnetic field replaces the spherical symmetry for an axisymmetrical one. 

As a first test, we consider the evolution of the star within the so-called "Cowling 
approximation", i.e. by holding the metric fixed to its initial value and by evolving the 
MHD variables onto this background spacetime (the evolution is not made only at the outer 
boundaries, where we use Dirichlet-type boundary conditions). The results of these evolutions 
are summarized in figure [8] with the left panel showing the evolution of the maximum 
of the rest-mass density p when normalized to its initial value. The three different lines 
(dotted, dashed and continuous) refer to the three resolutions used of N = 60 3 , 90 3 , and 
120 3 , respectively. The coordinate time on the horizontal axis is expressed in terms of the 
characteristic "dynamical timescale" r = J R? jM, where R is the coordinate radius of the 
star. 

In analogy with what observed in the purely hydrodynamical case 071 . the magnetized 
star is set into oscillation by the small truncation error introduced by the mapping 
onto a Cartesian coordinate system of the stationary solution found in spherical polar 
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Figure 8. Left panel: Maximum of the rest-mass density p normalized to its initial value and 
expressed in terms of the dynamical timescale r = *J Fft/M. The magnetized star is evolved 
within the Cowling approximation, with different lines referring to different resolutions: 
TV = 60 3 (dotted), TV = 90 3 (dashed line) and TV = 120 3 (solid line). Right panel: The 
same as in the left one but for a longer timescale. 



coordinates [73 1 (no perturbation coming from the outer boundaries was seen to influence the 
dynamics of the oscillations). Because of its stochastic nature, the initial perturbation triggers 
a variety of modes, most of which however decay rather rapidly leaving, after t ~ 25t, 
an essentially harmonic oscillation in the fundamental mode only. This is shown in the right 
panel of figure[8]which shows the evolution over a longer timescale. Furthermore, in the linear 
regime considered here, the amplitude of the oscillations is proportional to the magnitude of 
the truncation error and one expects the former to decrease as the resolution is increased. This 
is clearly the case for the oscillations shown in figure [8] and the very good overlap among the 
different timeseries is an indication that the oscillations indeed correspond to eigenmodes and 
that the code is tracking them correctly at these resolutions. 

Next, we consider the evolution of the same initial model discussed above but including 
also the solution of the Einstein equations so as to make the system fully dynamical (Dirichlet- 
type boundary conditions are used at the outer boundaries for the MHD variables and radiative 
ones for the fields). Also in this case, oscillations are triggered by the truncation error, with 
an amplitude that converges to zero with the increase of the resolution and with an harmonic 
content that becomes more evident after the initial transient (also in this case, no perturbation 
coming from the outer boundaries was seen to influence the dynamics of the oscillations). In 
addition, and in analogy with what observed in the purely hydrodynamical case Il74ll75ll37l . 
the oscillations are accompanied by a secular growth which also converges away at the correct 
rate with increasing grid resolution and that does not influence the long-term evolutions. This 
is shown in figure [9] which reports the same quantities as in figure|8]but for a fully dynamical 
evolution. Note also that the secular evolution of the central rest-mass density varies according 
to the EOS adopted: when using the ideal-fluid EOS, in fact, the secular drift of the central 
rest-mass density is towards lower densities. However, if the adiabatic condition is enforced, 
the opposite is true and central rest-mass density evolves towards larger values. Both the 
evolutions in the Cowling approximation and in dynamical spacetimes, have not shown signs 




Figure 9. Left panel: Maximum of the rest-mass density p normalized to its initial value and 
expressed in terms of the dynamical timescale r = sj i?, 3 /M. The magnetized star is evolved 
together with the spacetime, with different lines referring to different resolutions: N = 60 3 
(dotted), N = 90 3 (dashed line) and N = 120 3 (solid line). Right panel: The same as in the 
left one but for a longer timescale. 



of instability at all resolutions considered and up to several tens of dynamical timescales. 

As a final remark we underline that the convergence rate is not exactly second-order 
but slightly smaller, (i.e., 1.7-1.8), because the reconstruction schemes are only first-order 
accurate at local extrema {i.e. at the centre and at the surface of the star) thus increasing the 
overall truncation error. Similar estimates were obtained also in the purely hydrodynamical 
case |37|. 

5. Conclusions 

We have presented a new three-dimensional numerical code in Cartesian coordinates 
developed to solve the full set of GRMHD equation on a dynamical background. The 
code is based on high-resolution shock-capturing techniques as implemented on domains 
with adaptive mesh refinements. This code represents the extension to MHD of the 
approach already implemented with success in the general-relativistic hydrodynamics code 
Whisky l30l . 

The code has been validated through an extensive series of testbeds both in special and in 
general relativity scenarios. In particular, we have first considered a set of Riemann problems 
in a Minkowski spacetime following a variety of initial conditions. In all of the tests presented, 
the numerical solution has been compared with the exact one [49|, providing, for the first 
time, a quantitative assessment of the code's ability to evolve correctly and accurately all 
of the different waves that can form in relativistic MHD. Furthermore, as a demonstration 
of the proper handling of continuous and discontinuous flows in the presence of an excision 
region, we have extended the Riemann-problems tests across an excised boundary. In doing 
so we have revealed the importance of correct boundary conditions and pointed out that those 
recently proposed in [25 1 can lead to incorrect solutions for non-smooth flows. 

Next, to investigate magnetized fluids in a curved but fixed spacetime, we have 
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considered the spherical accretion of a perfect fluid with a radial magnetic field onto a 
Schwarzschild black hole (relativistic Bondi flow). Also in this case, the code has been shown 
to accurately reproduce the stationary solution and to be convergent at the correct rate for 
small and large magnetizations. For very large magnetizations, however, the very rapidly 
varying behaviour of the MHD variables near the excision region prevents from a correct 
convergence near the horizon, although the code remains second-convergent away from the 
horizon and is convergent overall. Also for these extremely-high values of the magnetic field, 
the code has shown to be robust and accurate at regimes where other codes were reported to 
fail 1251 . 

Finally, we have considered the evolution of magnetized neutron stars in equilibrium and 
constructed as a consistent solutions of the coupled Einstein-Maxwell equations. Such initial 
models represent an important difference from those considered by other authors, which were 
not consistent solution of the Einstein equations initially and whose magnetic field is totally 
confined within the star [25|. In analogy with what observed in the purely hydrodynamical 
case [37], these magnetized stars are set into oscillation by the small truncation error. These 
pulsations, which have been studied both in fixed (Cowling approximation) and in dynamical 
spacetimes, have been shown to have the correct behaviour under changes of spatial resolution 
and to correspond to the eigenmodes of relativistic and magnetized stars. Both evolutions in 
the Cowling approximation and in dynamical spacetimes have not shown signs of instability 
at all resolutions considered and up to several tens of dynamical timescales. 

A number of projects are expected to be carried out with the new code. Firstly, we plan 
to extend the study on the oscillations of rotating and nonrotating neutron stars with a detailed 
analysis of the effect of magnetic fields on the frequency of oscillations. It is important to note 
that only recently some results were obtained in perturbation theory and within the Cowling 
approximation [76|. Our code will be a complementary tool to the perturbative approaches, 
using the latter as testbeds and carrying them beyond the regimes of slow rotation and weak 
magnetizations. Such a study, and the comparison with the frequencies observed in objects 
such as the soft gamma-repeaters, will provide useful information on the mass and magnetic- 
field strength of magnetars. 

Secondly, we plan to use WhiskyMHD to study the collapse of uniformly and 
differentially rotating magnetized neutron stars with the aim of extending further the work 
done in |f37]|773 60 1 and to highlight the role that this process may have in the phenomenology 
of short 7-ray bursts. We are especially interested in the calculation of the gravitational-wave 
signal emitted by these sources and on the influence that magnetic fields may have on it. 

Finally, and in view of the total generality with which it has been developed, the code will 
be used to study the dynamics of binary systems of neutron stars and mixed binaries, with the 
aim of extending the work carried out in |78 1 and of considering in a fully general-relativistic 
context the Newtonian results obtained in 11241 . 

Acknowledgments 

It is a pleasure to thank L. Baiotti, G. Bodo, L. Del Zanna, I. Hawke, F. Loffler, C. D. Ott, E. 
Schnetter and O. Zanotti for useful discussions and comments. We are also indebted with J. 
Novak for his help in the construction of the equilibrium magnetized stars. All the numerical 
computations were performed on clusters Albert2 at the Physics Department of the University 
of Parma (Parma, Italy), CLX at CINECA (Bologna, Italy) and Peyote at the AEI (Golm, 
Germany). 



Whi skyMHD: a new numerical code for general relativistic magnetohydrodynamics 23 
References 

[1] J. R. Wilson. Some magnetic effects in stellar collapse and accretion. Ann. New York Acad. Sci., 262:123, 
1975. 

[2] J. A. Font. Numerical hydrodynamics in general relativity. Living Rev. Relativity, 6:4, 2003. URL (cited on 

December 15th 2006): http://www.livingreviews.org/lrr-2003-4. 
[3] M. Yokosawa. Energy and angular momentum transport in magnetohydrodynamical accretion onto a rotating 

black hole. Publ. Astron. Soc. Japan, 45:207-218, 1993. 
[4] M. Yokosawa. Structure and Dynamics of an Accretion Disk around a Black Hole. Publ. Astron. Soc. Japan, 

47:605-615, 1995. 

[5] S. Koide, K. Shibata, and T. Kudoh. Relativistic jet formation from black hole magnetized accretion disks: 
Method, tests, and applications of a general relativisticmagnetohydrodynamic numerical code. Astrophys 
J, 522:727-752, 1999. 

[6] S. F. Davis. A simplified tvd finite difference scheme via artificial viscosity. SIAMJ. Sci. Stat. Comput., 8:1-18, 
1987. 

[7] S. Koide, D. L. Meier, Shibata K., and Kudoh T. General relativistic simulations of early jet formation in a 

rapidly rotating black hole magnetosphere. Astrophys J., 536:668-674, 2000. 
[8] S. Koide, K. Shibata, T. Kudoh, and D. L. Meier. Extraction of black hole rotational energy by a magnetic field 

and the formation of relativistic jets. Science, 295:1688-1691, 2002. 
[9] S. Koide. Magnetic extraction of black hole rotational energy: Method and results of general relativistic 

magnetohydrodynamic simulations in kerr space-time. Phys. Rev. D, 67:104010, 2003. 
[10] J.-R De Villiers and J. F. Hawley. A numerical method for general relativistic magnetohydrodynamics. 

Astrophys. J, 589:458-480, 2003. 
[11] J. F. Hawley, J. R. Wilson, and L. L. Smarr. A numerical study of nonspherical black hole accretion, i equations 

and test problems. Astrophys. J., 277:296-311, 1984. 
[12] J. F. Hawley, J. R. Wilson, and L. L. Smarr. A numerical study of nonspherical black hole accretion, ii - finite 

differencing and code calibration. Astrophys. J., 55:211-246, 1984. 
[13] J.-R De Villiers, J. F. Hawley, and J. H. Krolik. Magnetically driven accretion flows in the kerr metric, i. models 

and overall structure. Astrophys. J, 599:1238-1253, 2003. 
[14] S. Hirose, J. H. Krolik, J.-R De Villiers, and J. F. Hawley. Magnetically driven accretion flows in the kerr 

metric, ii. structure of the magnetic field. Astrophys. J., 606:1083-1097, 2004. 
[15] J.-R De Villiers, J. F. Hawley, J. H. Krolik, and S. Hirose. Magnetically driven accretion in the kerr metric, iii. 

unbound outflows. Astrophys. J, 620:878-888, 2005. 
[16] C. F. Gammie, J. C. McKinney, and G. Toth. Harm: A numerical scheme for general relativistic 

magnetohydrodynamics. Astrophys. J, 589:444^157, 2003. 
[17] S. S. Komissarov. General relativistic magnetohydrodynamic simulations of monopole magnetospheres of 

black holes. MNRAS, 350:1431-1436, 2004. 
[18] L. Anton, O. Zanotti, J. A. Miralles, J. M. Marti, J. M. Ibanez, J. A. Font, and J. A. Pons. Numerical 3+1 

general relativistic magnetohydrodynamics: A local characteristic approach. Astrophys. J., 637:296-312, 

2006. 

[19] D. Neilsen, E. W. Hirschmann, and R. S. Millward. Relativistic MHD and excision: formulation and initial 

tests. Classical and Quantum Gravity, 23:505, 2006. 
[20] R J. Montero, O. Zanotti, J. A. Font, and L. Rezzolla. Dynamics of magnetized relativistic tori oscillating 

around black holes. MNRAS, submitted, 2007. 
[21] R D. Lax and B. Wendroff. Systems of conservation laws. Comm. Pure Appl. Math., 13:217-237, 1960. 
[22] T. Y. Hou and R G. LeFloch. Why Nonconservative Schemes Cconverge to Wrong Solutions: Error Analysis. 

Math. ofComp., 62:497-530, 1994. 
[23] M. Obergaulinger, M. A. Aloy, and E. Muller. Axisymmetric simulations of magneto-rotational core collapse: 

dynamics and gravitational wave signal. A&A, 450:1107-1134, 2006. 
[24] D. J. Price and S. Rosswog. Producing Ultrastrong Magnetic Fields in Neutron Star Mergers. Science, 

312:719-722,2006. 

[25] M. D. Duez, Y. T. Liu, S. L. Shapiro, and B. C. Stephens. Relativistic magnetohydrodynamics in dynamical 

spacetimes: Numerical methods and tests. Phys. Rev. D, 72:024028, 2005. 
[26] M. Shibata and Y. Sekiguchi. Magnetohydrodynamics in full general relativity: Formulation and tests. Phys. 

Rev. D, 72:044014, 2005. 

[27] M. D. Duez, Y. T. Liu, S. L. Shapiro, M. Shibata, and B. C. Stephens. Collapse of Magnetized Hypermassive 
Neutron Stars in General Relativity. Physical Review Letters, 96:031 101, 2006. 

[28] M. Shibata, M. D. Duez, Y. T. Liu, S. L. Shapiro, and B. C. Stephens. Magnetized Hypermassive Neutron-Star 
Collapse: A Central Engine for Short Gamma-Ray Bursts. Physical Review Letters, 96:031 102, 2006. 

[29] M. D. Duez, Y. T. Liu, S. L. Shapiro, M. Shibata, and B. C. Stephens. Evolution of magnetized, differentially 
rotating neutron stars: Simulations in full general relativity. Phys. Rev. D, 73:104015, 2006. 



WhiskyMHD: a new numerical code for general relativistic magnetohydrodynamics 



24 



[30] L. Baiotti, I. Hawke, P. Montero, and L. Rezzolla. A new three-dimensional general-relativistic hydrodynamics 
code. In R. Capuzzo-Dolcetta, editor, Computational Astrophysics in Italy: Methods and Tools, volume 1 , 
page 327, Trieste, 2003. Mem. Soc. Astron. It. Suppl. 

[31] Cactus, www.cactuscode.org. 

[32] T. Nakamura, K.-I. Oohara, and Y. Kojima. General relativistic collapse to black holes and gravitational waves 

from black holes. Prog. Theor. Phys. Suppl, 90:1-218, 1987. 
[33] M. Shibata and T. Nakamura. Evolution of 3-dimensional gravitational-waves - harmonic slicing case. Phys. 

Rev. D, 52:5428, 1995. 

[34] T. W. Baumgarte and S. L. Shapiro. Numerical integration of einstein's field equations. Phys. Rev. D, 
59:024007, 1999. 

[35] M. Alcubierre, B. Briigmann, T. Dramlitsch, JA. Font, P. Papadopoulos, E. Seidel, N. Stergioulas, and 
R. Takahashi. Towards a stable numerical evolution of strongly gravitating systems in general relativity: 
The conformal treatments. Phys. Rev. D, 62:044034, 2000. 

[36] M. Alcubierre, B. Briigmann, P. Diener, M. Koppitz, D. Pollney, E. Seidel, and R. Takahashi. Gauge conditions 
for long-term numerical black hole evolutions without excision. Phys. Rev. D, 67:084023, 2003. 

[37] L. Baiotti, I. Hawke, P. J. Montero, F. Loffler, L. Rezzolla, N. Stergioulas, J. A. Font, and E. Seidel. Three- 
dimensional relativistic simulations of rotating neutron-star collapse to a kerr black hole. Phys. Rev. D, 
71:024035, 2005. 

[38] J. M. Marti, J. M. Ibanez, and J. A. Miralles. Numerical relativistic hydrodynamics: local characteristic 

approach. Phys. Rev. D, 43:3794, 1991. 
[39] F. Banyuls, J. A. Font, J. M. Ibanez, J. M. Marti, and J. A. Miralles. Numerical 3+1 general relativistic 

hydrodynamics: A local characteristic approach. Astrophys. J., 476:221-231, 1997. 
[40] A. M. Anile. Relativistic Fluid and Magneto-Fluids. Cambridge University Press, 1989. 
[41] E. F. Toro. Riemann Solvers and Numerical Methods for Fluid Dynamics . Springer- Verlag, 1999. 
[42] S. A. Teukolsky. Stability of the iterated crank-nicholson method in numerical relativity. Phys. Rev. D, 

61:087501,2000. 

[43] G. Leiler and L. Rezzolla. Iterated crank-nicolson method for hyperbolic and parabolic equations in numerical 

relativity. Phys. Rev. D, 73:044001, 2006. 
[44] A. Harten, P. D. Lax, and B. van Leer. On upstream differencing and godunov-type schemes for hyperbolic 

conservation laws. SIAMRev, 25:35-61, 1983. 
[45] J. M. Marti and E. Muller. The analytical solution of the riemann problem in relativistic hydrodynamics. J. 

Fluid. Mech., 258:317, 1994. 
[46] J. A. Pons, J. M. Marti, and E. Muller. The exact solution of the riemann problem with non-zero tangential 

velocities in relativistic hydrodynamics. /. Fluid. Mech., 422:125, 2000. 
[47] L. Rezzolla and O. Zanotti. An improved exact riemann solver for relativistic hydrodynamics. /. Fluid. Mech., 

449:395, 2001. 

[48] L. Rezzolla, O. Zanotti, and J. A. Pons. An improved exact riemann solver for multi-dimensional relativistic 
flows. J. Fluid. Mech., 479:199, 2003. 

[49] B. Giacomazzo and L. Rezzolla. The exact solution of the riemann problem in relativistic magnetohydrody- 
namics. /. Fluid Mech., 562:223-259, 2006. 

[50] M. Abramowitz and A. Stegun. Handbook of Mathematical Functions . Dover, 1965. 

[51] W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, and M. Metcalf. Numerical Recipes in Fortran 

90. Cambridge University Press, 2nd edition, 2002. 
[52] T. Leismann, L. Anton, M. A. Aloy, E. Muller, J. M. Marti, J. A. Miralles, and J. M. Ibanez. Relativistic MHD 

simulations of extragalactic jets. A&A, 436:503-526, 2005. 
[53] J. U. Brackbill and D. C. Barnes. The effect of nonzero V ■ B on the numerical solution of the 

magnetohydrodynamic equations. /. Comp. Phys., 35:426^130, 1980. 
[54] K. S. Yee. Numerical solution of inital boundary value problems involving maxwell's equations in isotropic 

media. IEEE Trans. Antenna Propagation, AP- 14:302-307, 1966. 
[55] C. R. Evans and J. F. Hawley. Simulation of magnetohydrodynamic flows: A constrained transport method. 

Astrophys. J., 332:659, 1988. 
[56] D. S. Balsara and D. S. Spicer. A staggered mesh algorithm using high order godunov fluxes to ensure 

solenoidal magnetic fields in magnetohydrodynamic simulations. J. Comp. Phys., 149:270-292, 1999. 
[57] G. Toth. The V-6 = constraint in shock-capturing magnetohydrodynamics codes. /. Comp. Phys., 161:605- 

652, 2000. 

[58] S. C. Noble, C. F. Gammie, J. C. McKinney, and L. Del Zanna. Primitive variable solvers for conservative 
general relativistic magnetohydrodynamics. Astrophys. J., 641:626-637, 2006. 

[59] S. S. Komissarov. A godunov-type scheme for relativistic magnetohydrodynamics. MNRAS, 303:343-366, 
1999. 

[60] L. Baiotti and L. Rezzolla. Challenging the paradigm of singularity excision in gravitational collapse. Physical 
Review Utters, 97:141101, 2006. 



WhiskyMHD: a new numerical code for general relativistic magnetohydrodynamics 



25 



[61] E. Schnetter, S. H. Hawley, and I. Hawke. Evolutions in 3D numerical relativity using fixed mesh refinement. 

Class. Quantum Grav., 21(6):1465-1488, 2004. 
[62] M. J. Berger and J. Oliger. Adaptive mesh refinement for hyperbolic partial differential equations. /. Comput. 

Phys., 53:484-512, 1984. 
[63] Carpet, www.carpetcode.org. 

[64] D. S. Balsara. Total variation diminishing scheme for relativistic magnetohydrodynamics. ApJS, 132:83-101, 
2001. 

[65] M. Brio and C.-C. Wu. An upwind differencing scheme for the equations of ideal magnetohydrodynamics. J. 
Comp. Phys., 75:400, 1988. 

[66] M. H. P. M. Van Putten. A numerical implementation of mhd in divergence form. J. Comp. Phys., 105:339, 
1993. 

[67] D. Ryu and T. W. Jones. Numerical magnetohydrodynamics in astrophysics: Algorithm and tests for one- 
dimensional flow. Astrophys. /, 442:228, 1995. 

[68] S. A. E. G. Falle and S. S. Komissarov. On the inadmissibility of non-evolutionary shocks. J. Plasma Physics, 
65:29,2001. 

[69] M. Torrilhon. Non-uniform convergence of finite volume schemes for riemann problems of ideal 

magnetohydrodynamics. /. Comput. Phys., 192:73, 2003. 
[70] M. Torrilhon. Uniqueness conditions for riemann problems of ideal magnetohydrodynamics. J. Plasma 

Physics, 69:253, 2003. 

[71] M. Torrilhon and D. S. Balsara. High order weno schemes: investigations on non-uniform convergence for 

mhd riemann problems. /. Comput. Phys., 201:586, 2004. 
[72] F. C. Michel. Accretion of matter by condensed objects. Astrophys. Spa. Set, 15:153, 1972. 
[73] M. Bocquet, S. Bonazzola, E. Gourgoulhon, and J. Novak. Rotating neutron star models with a magnetic field. 

A&A, 301:757-775, 1995. 

[74] J. A. Font, N. Stergioulas, and K. D. Kokkotas. Nonlinear hydrodynamical evolution of rotating relativistic 

stars: Numerical methods and code tests. Mon. Not. R. Astron. Soc, 313:678, 2000. 
[75] J. A. Font, T. Goodale, S. Iyer, M. Miller, L. Rezzolla, E. Seidel, N. Stergioulas, W. M. Suen, and M. Tobias. 

Three-dimensional general relativistic hydrodynamics. II. Long-term dynamics of single relativistic stars. 

Phys. Rev. D, 65:084024, 2002. 
[76] H. Sotani, K. D. Kokkotas, and N. Stergioulas. Torsional Oscillations of Relativistic Stars with Dipole Magnetic 

Fields, astro-ph/0608626, 2006. 
[77] L. Baiotti, I. Hawke, L. Rezzolla, and E. Schnetter. Gravitational-Wave Emission from Rotating Gravitational 

Collapse in Three Dimensions. Physical Review Letters, 94:131 101, 2005. 
[78] F. LSffler, L. Rezzolla, and M. Ansorg. Numerical evolutions of a black hole-neutron star system in full General 

Relativity: Head-on collision. Phys Rev D, 74:104018, 2006. 



